################
## Preliminary Code
##################
## Set working directory equal to climate public opinion/DGIRT
setwd("/Users/Parrish/Dropbox (MIT)/Climate public opinion data/Replication")
library(dplyr)
library(lmtest)
library(foreign)
#library(plm)
library(tidyr)
library(lfe)
library(clusterSEs) ## for bootstrapped SEs


lag.new <- function(x, n = 1L, along_with){
 index <- match(along_with - n, along_with, incomparable = NA)
 out <- x[index]
 attributes(out) <- attributes(x)
 out
}

lead.new <- function(x, n = 1L, along_with){
 index <- match(along_with + n, along_with, incomparable = NA)
 out <- x[index]
 attributes(out) <- attributes(x)
 out
}

# Multiple plot function
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL) {
  require(grid)

                                        # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

                                        # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
                                        # Make the panel
                                        # ncol: Number of columns of plots
                                        # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), 
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }

  if (numPlots == 1) {
    print(plots[[1]])

  } else {
                                        # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

                                        # Make each plot, in the correct location
    for (i in 1:numPlots) {
                                        # Get the i, j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout ==  i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, 
                            layout.pos.col = matchidx$col))
    }
  }
}

clusterSE <- function(x, method="arellano", cluster="group") {
  if (any(grepl("plm", class(x)))) {
    coeftest(x, plm::vcovHC(x, method=method, cluster=cluster))
  }
}
plmHC <- function (mod) {
  plm::vcovHC(mod, method="arellano", cluster="group")
}
MOC <- function (data, its, model, vcov=plmHC, prog.int=100, rsq=TRUE) {
  ## Method of Composition (as described by Treier and Jackman)
  tildeB <- as.data.frame(matrix(nrow=length(its), ncol=length(coef(model))))
  R2 <- data.frame(rsq=rep(NA, length(its)), adjrsq=rep(NA, length(its)))
  for (s in seq_along(its)) {
    if (!s %% prog.int) print(s)
    it <- its[s]
    ## (1) Sample from p(x)
    data_s <- subset(data, It == it)
    names(data_s) <- gsub("\\_s$", "\\_mean", names(data_s)) 
    ## (2) Sample from p(B|x,y):
    ##     (a) Estimate B_s and Cov(B_s) conditional on x_s.
    mod_s <- update(model, data=data_s)
    hatB_s <- coef(mod_s)
    hatV_s <- vcov(mod_s)
    ##     (b) Sample \tilde{B_s} from MV(\hat{B_s}, \hat{Cov(B_s)}).
    tildeB[s, ] <- MASS::mvrnorm(n = 1, mu = hatB_s, Sigma = hatV_s)
    if (rsq) {
      R2$rsq[s] <- summary(mod_s)$r.squared["rsq"]
      R2$adjrsq[s] <- summary(mod_s)$r.squared["adjrsq"]
    }
  }
  names(tildeB) <- names(coef(model))
  if (rsq) {
    tildeB$rsq <- R2$rsq
    tildeB$adjrsq <- R2$adjrsq
  }
  return(tildeB)
}
pEmp <- function (x) {
  p_pos <- mean(x < 0)
  p_neg <- mean(x > 0)
  p2side <- 2 * min(p_pos, p_neg)
  p2side
}
pNorm <- function (x) {
  z <- abs(mean(x)/sd(x))
  p2side <- 2*(1 - pnorm(z))
  p2side
}
MOCsumm <- function (moc, regex="tmpc_lagged|concern |tmax_lagged |treatment|Intercept|rsq|mass_climate_concern_lagged", digits=3,
                     ff=funs(est = mean, se = sd, z = mean(.)/sd(.),
                             pnorm = pNorm(.),
                             pemp = pEmp(.))) {
  summarise_at(moc, .vars=vars(matches(regex)), .funs=ff) %>%
    mutate_all(funs(round(., digits=digits))) %>%
    reshape2::melt(measure.vars=names(.))
}


MOCmean <- function (mod.ls) {
  if (all(plyr::laply(mod.ls, function (x) "moc" %in% names(x)))) {
    plyr::llply(mod.ls, function (x) apply(x$moc, 2, mean))
  } else {
    plyr::llply(mod.ls, function (x) apply(x, 2, mean))
  }
}
MOCsd <- function (mod.ls) {
  if (all(plyr::laply(mod.ls, function (x) "moc" %in% names(x)))) {
    plyr::llply(mod.ls, function (x) apply(x$moc, 2, sd))
  } else {
    plyr::llply(mod.ls, function (x) apply(x, 2, sd))
  }
}

FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}

mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}

## flag for whether to do robustness checks that aren't in main paper
extra.analysis<-0

### Output the main dataset 
climate_scale2<- read.dta(file="climate_analysis_states.dta")

climate_scale2$state_lower<-tolower(climate_scale2$NAME)

### standardize the various climate measures

climate_scale2<-subset(climate_scale2, year>1998)

climate_scale2$log_events_lagged<-log(climate_scale2$events_lagged+1)
climate_scale2$log_wildfire_lagged<-log(climate_scale2$wildfire_lagged+1)


climate_scale2$tmpc_lagged_std<-scale(climate_scale2$tmpc_lagged)[,1]
climate_scale2$tmax_lagged_std<-scale(climate_scale2$tmax_lagged)[,1]
climate_scale2$tmin_lagged_std<-scale(climate_scale2$tmin_lagged)[,1]
climate_scale2$hdd_lagged_std<-scale(log(climate_scale2$hdd_lagged))[,1]
climate_scale2$cdd_lagged_std<-scale(log(climate_scale2$cdd_lagged))[,1]
climate_scale2$log_events_lagged_std<-scale(climate_scale2$log_events)[,1]
climate_scale2$pdsi_lagged_std<-scale(climate_scale2$pdsi_lagged)[,1]
climate_scale2$phdi_lagged_std<-scale(climate_scale2$phdi_lagged)[,1]
climate_scale2$sp12_lagged_std<-scale(climate_scale2$sp12_lagged)[,1]
climate_scale2$log_wildfire_lagged_std<-scale(climate_scale2$log_wildfire_lagged)[,1]



### run vcov cluster ###
#library(lmtest)
climate_scale2<-subset(climate_scale2, year>1988)
climate_scale2$abb<-as.vector(climate_scale2$abb)

##########
## Main empirical analysis
##########

climate_scale3<-climate_scale2

climate_scale3<-dplyr::group_by(climate_scale3, year, abb) %>%
	dplyr::summarise_each(funs(mean) )
	
climate_scale3 <- dplyr::group_by(climate_scale3, abb) %>%
  dplyr::mutate(mass_climate_concern_lagged = lag.new(mass_climate_concern, 1, along_with = year)   )
  
  
climate_scale3$time<-climate_scale3$year-1998
climate_scale3$time_sq<-climate_scale3$time^2


##########
## Table F1: Non-Measurement Error Adjusted Main Empirical Analysis
##########

##########
## models with state fixed effect and within state trend
##########
## Column 1
mod1.temp <- felm(mass_climate_concern ~ tmpc_lagged+abb*time | abb|0|abb,data=climate_scale3)
summary(mod1.temp)

##########
## models with two-way fixed effect
##########
mod2.temp <- felm(mass_climate_concern ~ tmpc_lagged | abb+year|0|abb,data=climate_scale3)
summary(mod2.temp)

####
## Table G2: Models by time period, non-measurement error adjusted
####
mod2.temp.early <- felm(mass_climate_concern ~ tmpc_lagged | abb+year|0|abb,data=climate_scale3[climate_scale3$year<2005,])
summary(mod2.temp.early)

mod2.temp.mid <- felm(mass_climate_concern ~ tmpc_lagged | abb+year|0|abb,data=climate_scale3[climate_scale3$year>2004 & climate_scale3$year<2011,])
summary(mod2.temp.mid)

mod2.temp.late <- felm(mass_climate_concern ~ tmpc_lagged | abb+year|0|abb,data=climate_scale3[climate_scale3$year>2010,])
summary(mod2.temp.late)

library(stargazer)
stargazer(mod2.temp.early, mod2.temp.mid, mod2.temp.late)

##########
## models with fixed effects for state and linear time trends in each state
##########

## Column 3
mod3.temp <- felm(mass_climate_concern ~ tmpc_lagged+abb*year | abb+year|0|abb,data=climate_scale3)
summary(mod3.temp)


##########
## models with state fixed effects and LDV
##########
## Column 4
mod4.temp <- felm(mass_climate_concern ~ tmpc_lagged+mass_climate_concern_lagged | abb+year|0|abb,data=climate_scale3)
summary(mod4.temp)

library(stargazer)
stargazer(mod2.temp, mod3.temp,mod4.temp)

##########
## Robustness to different measures of temperature
##########

mod2a.temp <- felm(mass_climate_concern ~ tmpc_lagged_std | abb+year|0|abb,data=climate_scale3)
summary(mod2a.temp)

mod2b.temp <- felm(mass_climate_concern ~ tmax_lagged_std | abb+year|0|abb,data=climate_scale3)
summary(mod2b.temp)

mod2c.temp <- felm(mass_climate_concern ~ tmin_lagged_std | abb+year|0|abb,data=climate_scale3)
summary(mod2c.temp)

mod2d.temp <- felm(mass_climate_concern ~ hdd_lagged_std | abb+year|0|abb,data=climate_scale3)
summary(mod2d.temp)

mod2e.temp <- felm(mass_climate_concern ~ cdd_lagged_std | abb+year|0|abb,data=climate_scale3)
summary(mod2e.temp)
      
stargazer(mod2a.temp,mod2b.temp,mod2c.temp,mod2d.temp,mod2e.temp)
  
##########
## Graph of non-measurement error corrected regression results (not in paper)
##########

if(extra.analysis==1){

treatment<- c("tmpc_lagged_std","log_events_lagged_std", "pdsi_lagged_std", "phdi_lagged_std", "sp12_lagged_std", "log_wildfire_lagged_std")


library(dotwhisker)
library(broom)
library(dplyr)
mod1.temp <- NULL
mod2.temp <- NULL
mod3.temp <- NULL
mod4.temp <- NULL
mod5.temp <- NULL
mod6.temp <- NULL


mod1.temp.cl <- NULL
mod2.temp.cl <- NULL
mod3.temp.cl <- NULL
mod4.temp.cl <- NULL
mod5.temp.cl <- NULL
mod6.temp.cl <- NULL

for (p in 1:length(treatment)){
print(p)

##########
## models with state fixed effects 
##########

climate_scale3$year<-as.numeric(as.vector(climate_scale3$year))
## Column 1
mod1.temp[[p]] <- felm(mass_climate_concern ~ (get(treatment[p])) | abb|0|abb,data=climate_scale3)
summary(mod1.temp[[p]] )
mod1.temp.cl[[p]]<-tidy(summary(mod1.temp[[p]])[7]$coefficients)
colnames(mod1.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")

##########
## models with two-way fixed effects
##########

mod2.temp[[p]] <- felm(mass_climate_concern ~ (get(treatment[p])) | abb+year|0|abb,data=climate_scale3)
summary(mod2.temp[[p]] )

mod2.temp.cl[[p]]<-tidy(summary(mod2.temp[[p]])[7]$coefficients)
colnames(mod2.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")

##########
## models with fixed effects for state and linear time trends in each state
##########

## Column 3
mod3.temp[[p]] <- felm(mass_climate_concern ~ (get(treatment[p]))+abb*year | abb+year|0|abb,data=climate_scale3)
summary(mod3.temp[[p]] )
mod3.temp.cl[[p]]<-tidy(summary(mod3.temp[[p]])[7]$coefficients)
colnames(mod3.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")

##########
## models with state fixed effects and LDV
##########


mod4.temp[[p]] <- felm(mass_climate_concern ~ (get(treatment[p]))+mass_climate_concern_lagged | abb+year|0|abb,data=climate_scale3)
summary(mod4.temp[[p]] )
mod4.temp.cl[[p]]<-tidy(summary(mod4.temp[[p]])[7]$coefficients)
colnames(mod4.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")
}

results_unadjusted<-NULL
for (p in 1:length(treatment)){
	results_unadjusted<-rbind(results_unadjusted,mod2.temp.cl[[p]][1,],mod3.temp.cl[[p]][1,],mod4.temp.cl[[p]][1,])
}

results_unadjusted$model<-rep(c("(1) State/Year FEs", "(2): State/Year FEs + Time Trend", "(3) State/Year FEs + LDV"),6)

results_unadjusted$term<-c(rep("Average Temperature",3),rep("Storm Episodes",3),rep("Short-term Drought Severity Index",3),rep("Long-term Drought Severity Index",3),rep("Precipitation Index",3),rep("Wildfires",3))

n_vars2 <- length(unique(results_unadjusted$term))
n_models2 <- length(unique(results_unadjusted$model))
y_ind2<-rep(seq(n_vars2, 1), n_models2)
var_names2<-results_unadjusted$term
pdf('ExtraFigure_regression_results_unadjusted.pdf', height = 4.25)
dwplot(results_unadjusted, fill=model, alpha=.1) +
     theme_bw() + xlab("Coefficient Estimate") + ylab("") +
     geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
 #    ggtitle("Effect of Various Climate Indicators
  #    on Public Opinion")  +
     theme(plot.title = element_text(face="bold"), legend.position="bottom",
           legend.background = element_rect(colour="grey80"),legend.text = element_text(colour="black", size=7))+ theme(legend.title=element_blank())
   dev.off()  
}

##########
## Measurement Error Adjusted Empirical Analysis
##########
climate_scale3$time<-climate_scale3$year-1998
climate_scale3$time_sq<-climate_scale3$time^2

climate_scale4<-climate_scale2
#climate_scale4$It<-climate_scale4$iteration
#its.keep<-sample(1:1252, 400)
#climate_scale4<-climate_scale4[climate_scale4$It  %in% its.keep,]
climate_scale4$time<-climate_scale4$year-1998
climate_scale4$time_sq<-climate_scale4$time^2

climate_scale4 <- dplyr::group_by(climate_scale4, abb, It) %>%
  dplyr::mutate(mass_climate_concern_lagged = lag.new(mass_climate_concern, 1, along_with = year)   )

its<- unique(climate_scale4$It)
#its <- sort(sample(unique(climate_scale4$It), min(400, length(unique(climate_scale4$It)))))


### PANEL DATA FRAME
library(plm)
climate_scale.pd.summ <- pdata.frame(climate_scale3, index=c("abb", "year"), row.names=FALSE)
climate_scale.pd <- pdata.frame(climate_scale4, index=c("abb", "year"), row.names=FALSE)


model.names <- c(
  ## Table 1
 "oneway", "fe", "fe.tt","fe.ldv" )
model.ls <- vector("list", length=length(model.names))
names(model.ls) <- model.names

## one-way FE, measurement error adjusted
model.ls[["oneway"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb*time+abb,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["oneway"]]$coeftest <-
   clusterSE(model.ls[["oneway"]]$model))
model.ls[["oneway"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["oneway"]]$model)
MOCsumm(model.ls[["oneway"]]$moc)
model.ls[["oneway"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
  
## Table 1, column 1
## two-way FE, measurement error adjusted
model.ls[["fe"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe"]]$coeftest <-
   clusterSE(model.ls[["fe"]]$model))
model.ls[["fe"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe"]]$model)
MOCsumm(model.ls[["fe"]]$moc)
model.ls[["fe"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
  
## Table 1, column 2
## two-way FE with state time trends, measurement error adjusted
model.ls[["fe.tt"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb*time+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.tt"]]$coeftest <-
   clusterSE(model.ls[["fe.tt"]]$model))
model.ls[["fe.tt"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.tt"]]$model)
MOCsumm(model.ls[["fe.tt"]]$moc)
model.ls[["fe.tt"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
 
## Table 1, column 3
model.ls[["fe.ldv"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+mass_climate_concern_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.ldv"]]$coeftest <-
   clusterSE(model.ls[["fe.ldv"]]$model))
model.ls[["fe.ldv"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.ldv"]]$model)
MOCsumm(model.ls[["fe.ldv"]]$moc)
model.ls[["fe.ldv"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
    

## Table G1, column 1
## two-way FE, measurement error adjusted
model.ls[["fe.early"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb+year,
      data=climate_scale.pd.summ[as.numeric(as.vector(climate_scale.pd.summ$year))<2005,], model="pooling")
(model.ls[["fe.early"]]$coeftest <-
   clusterSE(model.ls[["fe.early"]]$model))
model.ls[["fe.early"]]$moc <-
  MOC(data=climate_scale.pd[as.numeric(as.vector(climate_scale.pd$year))<2005,], its=its, model.ls[["fe.early"]]$model)
MOCsumm(model.ls[["fe.early"]]$moc)
model.ls[["fe.early"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
  
## Table G1, column 2
## two-way FE, measurement error adjusted
model.ls[["fe.mid"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb+year,
      data=climate_scale.pd.summ[as.numeric(as.vector(climate_scale.pd.summ$year))>2004 &as.numeric(as.vector(climate_scale.pd.summ$year))<2011,], model="pooling")
(model.ls[["fe.mid"]]$coeftest <-
   clusterSE(model.ls[["fe.mid"]]$model))
model.ls[["fe.mid"]]$moc <-
  MOC(data=climate_scale.pd[as.numeric(as.vector(climate_scale.pd$year))>2004 &as.numeric(as.vector(climate_scale.pd$year))<2011,], its=its, model.ls[["fe.mid"]]$model)
MOCsumm(model.ls[["fe.mid"]]$moc)
model.ls[["fe.mid"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()

## Table G1, column 3
## two-way FE, measurement error adjusted
model.ls[["fe.late"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb+year,
      data=climate_scale.pd.summ[as.numeric(as.vector(climate_scale.pd.summ$year))>2010,], model="pooling")
(model.ls[["fe.late"]]$coeftest <-
   clusterSE(model.ls[["fe.late"]]$model))
model.ls[["fe.late"]]$moc <-
  MOC(data=climate_scale.pd[as.numeric(as.vector(climate_scale.pd$year))>2010,], its=its, model.ls[["fe.late"]]$model)
MOCsumm(model.ls[["fe.late"]]$moc)
model.ls[["fe.late"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
  

if(extra.analysis==1){

#############
## Robustness checks: Not in paper 
##########

model.ls[["fe.average"]]$model <-
 plm(mass_climate_concern ~
        tmpc_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.average"]]$coeftest <-
   clusterSE(model.ls[["fe.average"]]$model))
model.ls[["fe.average"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.average"]]$model)
MOCsumm(model.ls[["fe.average"]]$moc)
model.ls[["fe.average"]]$moc %>%
  ggplot(aes(x=tmpc_lagged)) +
  geom_histogram()
    
model.ls[["fe.max"]]$model <-
 plm(mass_climate_concern ~
        tmax_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.max"]]$coeftest <-
   clusterSE(model.ls[["fe.max"]]$model))
model.ls[["fe.max"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.max"]]$model)
MOCsumm(model.ls[["fe.max"]]$moc)
model.ls[["fe.max"]]$moc %>%
  ggplot(aes(x=tmax_lagged)) +
  geom_histogram()

model.ls[["fe.min"]]$model <-
 plm(mass_climate_concern ~
        tmin_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.min"]]$coeftest <-
   clusterSE(model.ls[["fe.min"]]$model))
model.ls[["fe.min"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.min"]]$model)
MOCsumm(model.ls[["fe.min"]]$moc)
model.ls[["fe.min"]]$moc %>%
  ggplot(aes(x=tmin_lagged)) +
  geom_histogram()

}
### ##########
## Figure 4: Graph of measurement error corrected regression results
##########

treatment<- c("tmpc_lagged_std","log_events_lagged_std", "pdsi_lagged_std", "phdi_lagged_std", "sp12_lagged_std", "log_wildfire_lagged_std")

models <- lapply(treatment, function(x) {

   lm(substitute(mass_climate_concern ~ i+time+factor(abb), list(i = as.name(x))), data = climate_scale3)
})

lapply(models, summary)[1]

cor(climate_scale3$tmpc_lagged_std, climate_scale3$pdsi_lagged_std)

library(dotwhisker)
library(broom)
library(dplyr)
mod1.temp <- NULL
mod2.temp <- NULL
mod3.temp <- NULL
mod4.temp <- NULL
mod5.temp <- NULL
mod6.temp <- NULL


mod1.temp.cl <- NULL
mod2.temp.cl <- NULL
mod3.temp.cl <- NULL
mod4.temp.cl <- NULL
mod5.temp.cl <- NULL
mod6.temp.cl <- NULL


for (p in 1:length(treatment)){
print(p)


model.ls[["fe"]]$model <-
 plm(mass_climate_concern ~
       (get(treatment[p]))+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe"]]$coeftest <-
   clusterSE(model.ls[["fe"]]$model))
model.ls[["fe"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe"]]$model)
out<-MOCsumm(model.ls[["fe"]]$moc)

mod2.temp.cl[[p]]<-tidy(cbind(treatment[p], out[2,2], out[6,2] , out[10,2], out[14,2]))
colnames(mod2.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")

model.ls[["fe.tt"]]$model <-
 plm(mass_climate_concern ~
       (get(treatment[p]))+abb*time+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.tt"]]$coeftest <-
   clusterSE(model.ls[["fe.tt"]]$model))
model.ls[["fe.tt"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.tt"]]$model)
out<-MOCsumm(model.ls[["fe.tt"]]$moc)

mod3.temp.cl[[p]]<-tidy(cbind(treatment[p], out[2,2], out[6,2] , out[10,2], out[14,2]))
colnames(mod3.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")
    

model.ls[["fe.ldv"]]$model <-
 plm(mass_climate_concern ~
       (get(treatment[p]))+mass_climate_concern_lagged+abb+year,
      data=climate_scale.pd.summ, model="pooling")
(model.ls[["fe.ldv"]]$coeftest <-
   clusterSE(model.ls[["fe.ldv"]]$model))
model.ls[["fe.ldv"]]$moc <-
  MOC(data=climate_scale.pd, its=its, model.ls[["fe.ldv"]]$model)
out<-MOCsumm(model.ls[["fe.ldv"]]$moc)

mod4.temp.cl[[p]]<-tidy(cbind(treatment[p], out[2,2], out[7,2] , out[12,2], out[17,2]))
colnames(mod4.temp.cl[[p]])<-c("term", "estimates", "std.error", "statistic", "p.value")
    
}

results<-NULL
for (p in 1:length(treatment)){
	results<-rbind(results,mod2.temp.cl[[p]][1,],mod3.temp.cl[[p]][1,],mod4.temp.cl[[p]][1,])
#	results<-rbind(results,mod2.temp.cl[[p]][2,])
}

results$model<-rep(c("(1) State/Year FEs", "(2) State/Year FEs + Time Trend", "(3) State/Year FEs + LDV"),6)
#results$model<-rep(c("(1) State/Year FEs"),6)

results$term<-c(rep("Average Temperature",3),rep("Storm Episodes",3),rep("Short-term Drought Severity Index",3),rep("Long-term Drought Severity Index",3),rep("Precipitation Index",3),rep("Wildfires",3))

n_vars2 <- length(unique(results$term))
n_models2 <- length(unique(results$model))
y_ind2<-rep(seq(n_vars2, 1), n_models2)
var_names2<-results$term
pdf('Figure4_regression_results_measurement_error_adjusted.pdf', height = 4.25)
dwplot(results, fill=model, alpha = .1) +
     theme_bw() + xlab("Coefficient Estimate") + ylab("") +
     geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
 #    ggtitle("Effect of Various Climate Indicators
  #    on Public Opinion")  +
     theme(plot.title = element_text(face="bold"), legend.position="bottom",
           legend.background = element_rect(colour="grey80"),legend.text = element_text(colour="black", size=7))+ theme(legend.title=element_blank())
   dev.off()  

#####
### Figure C2: Graph national temperature
######   


### Non-measurement error adjusted analysis at national level
national_temperature2<- read.dta(file="climate_analysis_national.dta")

national_temperature3<-national_temperature2

national_temperature3<-dplyr::group_by(national_temperature2, year) %>%
	dplyr::summarise_each(funs(mean) )
	
national_temperature_quantiles<-dplyr::group_by(national_temperature2, year) %>%
	summarise(
		mass_climate_concern_lower=quantile(mass_climate_concern, probs=0.05),
		mass_climate_concern_median=quantile(mass_climate_concern, probs=0.5),
		mass_climate_concern_upper=quantile(mass_climate_concern, probs=0.95)
	)
national_temperature3<-merge(national_temperature3, national_temperature_quantiles, by="year")	
gallup_worry<-read.csv(file="gallup_summary.csv")
climate_scale_national<-merge(national_temperature3, gallup_worry, by.x="year",by.y="yr", all.x=T)
cor(climate_scale_national$mass_climate_concern_median, climate_scale_national$worry, use="complete.obs")
cor(climate_scale_national$mass_climate_concern_median, climate_scale_national$average_temperature, use="complete.obs")
climate_scale_national$time<-climate_scale_national$year-1998


 ######   
 r<-cor(climate_scale_national$mass_climate_concern_median, climate_scale_national$worry, use="complete.obs")
p <- ggplot(climate_scale_national, aes(x = mass_climate_concern_median, y = worry, label=year)) 
p <- p + geom_errorbar(aes(ymin=worry-1.96*worry_se,ymax=worry+1.96*worry_se,width=0.01), colour="grey60")
p <- p + geom_text()
p <- p + 
     theme_bw()+ ylab("Percent Worried about Climate Change on Gallup Poll")+ xlab("Climate Opinion Index")
p <- p + stat_smooth(method = "lm", formula = y ~ x, size = 1)
p <- p + annotate("text",x=-.9, y=.7, label=paste("r=", round(r,2), sep=""), size=9)
pdf('FigureC2_national_opinion_gallup_worry.pdf', height = 6)
print(p)
dev.off()
   

national<-(lm(mass_climate_concern_median~tmpc_lagged, data=climate_scale_national))
national_trends<-(lm(mass_climate_concern_median~tmpc_lagged+time, data=climate_scale_national))
stargazer(national,national_trends)

### Table E1: Measurement error adjusted analysis at national level

its<- unique(national_temperature2$It)

climate_scale.national.summ <- pdata.frame(national_temperature3, index=c("year"), row.names=FALSE)
climate_scale.national.pd <- pdata.frame(national_temperature2, index=c("year"), row.names=FALSE)

model.ls[["national"]]$model <-
 plm(mass_climate_concern ~
      tmpc_lagged,
      data=climate_scale.national.summ, model="pooling")
model.ls[["national"]]$moc <-
  MOC(data=climate_scale.national.pd, its=its, model.ls[["national"]]$model)
out<-MOCsumm(model.ls[["national"]]$moc)

 
 
#####
### Figure 2 in the main paper
#####

a <- ggplot(climate_scale_national, aes(x = year, y = mass_climate_concern_median)) 
a <- a + geom_ribbon(aes(ymin=mass_climate_concern_lower,ymax=mass_climate_concern_upper,width=0.2), alpha=0.2)
a <- a + geom_point()
a <- a + scale_y_continuous(breaks = c(-2,-1,0,1,2,3), labels=c("   -2","   -1","    0","    1","    2","    3"))
a <- a + scale_x_continuous(breaks = c(2000,2005,2010, 2015), limits=c(1998, 2017))
a <- a + geom_line()+
     theme_bw()+ ylab("Climate Concern Index")+ xlab("Year")

b<-(ggplot(climate_scale_national, aes(x=year, y=average_temperature))
  + geom_line(data=subset(climate_scale_national, year>1997 &  year<2018),aes(x=year, y=average_temperature),se = FALSE)
  + geom_point(data=subset(climate_scale_national, year>1997 &  year<2018),aes(x=year, y=average_temperature), size = rel(2))
  + scale_x_continuous(breaks = c(2000,2005,2010, 2015), limits=c(1998, 2017))
  + labs(x = 'Year', y = 'National Ave. Temp. (C)')
  + theme_bw() 
  + theme(legend.position = 'bottom')
 )
 b

pdf('Figure2_temp_opinion.pdf', height = 4)
multiplot(a, b)
dev.off()
   
save.image("climate_results.RData")